# -*- coding: utf-8 -*-
"""
Created on Thu Jul 28 08:18:36 2022

@author: Derek Rodriguez Pacheco, PhD Candidate
IUSS Pavia
"""

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import subprocess
import matplotlib as mpl
from Spectra.SpectrumA import spectrumAD 

mpl.rc('font',family = 'Times New Roman')

m= 960                           ## Mass in N*s2/m or kg
g= 9.81                          ## Acceleration of gravity in m/s2
f1 = 3.25
k1 = (4*np.pi**2)*f1**2*m

adjfactors= np.array(pd.read_csv('Adjustment_Factors_Displacement_CM.txt', header= None, delim_whitespace= True)[0])
indexesb= np.array(pd.read_csv('Indexes_Beginning.txt', header= None, delim_whitespace= True)[0])
indexesf= np.array(pd.read_csv('Indexes_End.txt', header= None, delim_whitespace= True)[0])

## Model's parameters

k1= 0.40
k2= 0.25
sigAct= 0.6
betaf= 1.3
epsSlip= 0
epsBear= 0
rBear= 1

testnumber= [13,14,15,16,17,18,19]

##

datacm= {}
# acceltcm01= {}
# acceltcm02= {}
disp1 = {}
disp2 = {}
dispCM = {}
accelcm= {}
timecm= {}
timecmc= {}

data= {}
datad= {}
accelda= {}
# accelda1= {}
# accelro= {}
accelbase= {}
dispda= {}
dispdacm= {}
timearray= {}
energyarray= {}
energyarraycm= {}
maxaccel= []
minaccel= []
maxaccelf= []
matchindexmin= []
matchindexmax= []
matchdispmin= []
matchdispmax= []
timetotal= []
acceltotal= []
dispnet= {}
dispnetcm= {}
disp1 = {}
disp2 = {}
# timestart= []

for i in range(len(testnumber)):
    datacm['%d' % (i+1)]= np.array(pd.read_csv('Data_Center_Mass\\accelAndDispl_15Hz\\%d_AccelAndDispl_15Hz.txt' % testnumber[i], header= None, delim_whitespace= True))
    accelcm['%d' % (i+1)]= np.array(datacm['%d' % (i+1)][:,4])[indexesb[i]:indexesf[i]+1]*(-1)
    timecm['%d' % (i+1)]= np.array(datacm['%d' % (i+1)][:,0])[indexesb[i]:indexesf[i]+1]
    timecmc['%d' % (i+1)]= timecm['%d' % (i+1)]-timecm['%d' % (i+1)][0]
    disp1['%d' % (i+1)] = np.array(datacm['%d' % (i+1)][:,5])[indexesb[i]:indexesf[i]+1]
    disp2['%d' % (i+1)] = np.array(datacm['%d' % (i+1)][:,6])[indexesb[i]:indexesf[i]+1]

for i in range(len(testnumber)):
    dm = []
 
    for j in range(len(disp1['%d' % (i+1)])):
        d1 = disp1['%d' % (i+1)][j]
        d2 = disp2['%d' % (i+1)][j]
        dm.append((d1+d2)/2)
        
    dispdacm['%d' % (i+1)] = np.array(dm)*adjfactors[i]

data['1']= np.array(pd.read_csv('13_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['2']= np.array(pd.read_csv('14_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['3']= np.array(pd.read_csv('15_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['4']= np.array(pd.read_csv('16_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['5']= np.array(pd.read_csv('17_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['6']= np.array(pd.read_csv('18_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['7']= np.array(pd.read_csv('19_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))



for i in range(len(testnumber)):
    accelbase['%d' % (i+1)]= np.array(data['%d' % (i+1)][:,1])
    timearray['%d' % (i+1)]= data['%d' % (i+1)][:,0]-data['%d' % (i+1)][0,0]
    energyarraycm['%d' % (i+1)]= [np.trapz(accelcm['%d' % (i+1)][:j]*m,dispdacm['%d' % (i+1)][:j]/1000) for j in range(len(dispdacm['%d' % (i+1)]))]
    
cc= 0

for i in range(len(testnumber)):
    for j in range(len(timearray['%d' % (i+1)])):
        if i==0:
            timetotal.append(timearray['%d' % (i+1)][j])
        else:
            timetotal.append(cc+timearray['%d' % (i+1)][j])
        acceltotal.append(accelbase['%d' % (i+1)][j])
    cc= timetotal[-1]
    if i!=len(testnumber)-1:
        for j in range(int(round(15/timearray['%d' % (i+1)][1]))):
            if i==0:
                timetotal.append(timearray['%d' % (i+1)][-1]+(j+1)*timearray['%d' % (i+1)][1])
            else:
                timetotal.append(cc+(j+1)*timearray['%d' % (i+1)][1])
            acceltotal.append(0)
        cc= timetotal[-1]

#### SDOF Analysis

# for i in range(15):
np.savetxt('TH_SDOF/Time_Total.txt', timetotal)
np.savetxt('TH_SDOF/Acceleration_Total.txt', acceltotal)

recordsa= 'TH_SDOF/Acceleration_Total.txt'
recordst= 'TH_SDOF/Time_Total.txt'

dts= 0.001953
mass= m/1e6

dur= timetotal[-1]

fh= open('FMdataSDOF.tcl', 'w')

fh.write(
     'set FMfilet "%s";\nset FMfilea "%s";\nset dtg %f;\nset durs %f;'
    '\nset M %f;'
    '\nset k1 %f;\nset k2 %f;\nset sigAct %f;\nset betaf %f;\nset epsSlip %f;\nset epsBear %f;\nset rBear %f;' % (recordst, recordsa, dts, dur, mass, k1, k2, sigAct, betaf, epsSlip, epsBear, rBear))
    # % (recordst, recordsa, dts, dur, mass, s1p, e1p, s2p, e2p, s3p, e3p, pinchX, pinchY, damage1, damage2, beta, s1n, e1n, s2n, e2n, s3n, e3n, width, height))
fh.close()
subprocess.call('OpenSees SDOF_CalibrationMH_Grav.tcl', shell=True)

####

sdofd= {}
sdofa= {}
sdoft= {}
sdoff= {}
sdofe= {}

dtt = 0.001953

indexes_Beg = [0,int(42/dtt),int(85/dtt),int(127/dtt),int(170/dtt),int(212/dtt),int(255/dtt)]
indexes_End = [int(28/dtt),int(71/dtt),int(113/dtt),int(156/dtt),int(198/dtt),int(241/dtt),144914]
timei= [0,42,85,127,170,212,255]

for i in range(len(testnumber)):
    sdofd['%d' % (i+1)]= np.array(pd.read_csv('ResultsAttSDOF/DispT.out', header= None, delim_whitespace= True)[1][indexes_Beg[i]:indexes_End[i]])*0.001
    sdofa['%d' % (i+1)]= np.array(pd.read_csv('ResultsAttSDOF/AccelT.out', header= None, delim_whitespace= True)[1][indexes_Beg[i]:indexes_End[i]])
    sdoft['%d' % (i+1)]= np.array(pd.read_csv('ResultsAttSDOF/AccelT.out', header= None, delim_whitespace= True)[0][indexes_Beg[i]:indexes_End[i]])-timei[i]
    sdoff['%d' % (i+1)]= sdofa['%d' % (i+1)]*0.001*m*-1
    sdofe['%d' % (i+1)]= [np.trapz(sdoff['%d' % (i+1)][:j],sdofd['%d' % (i+1)][:j]) for j in range(len(sdoff['%d' % (i+1)]))]


## Plotting

nr= 3
nc= 3

fig1, ax1 = plt.subplots(nr,nc,figsize=(13,8))

for i in range(nr-1):
    for j in range(nc):
        ns = int(3*i+j+1)
        ax1[i,j].set_facecolor('#DEE2E6')
        ax1[i,j].grid(color= 'black', ls= '--', lw= 0.5)
        ax1[i,j].set_xticks([0,10,20,30])
        ax1[i,j].set_yticks([0,100,200,300,400,500,600,700,800])
        ax1[i,j].set_xlim([0,30])
        ax1[i,j].set_ylim([0,800])
        ax1[i,j].plot(timecmc['%d' % ns],energyarraycm['%d' % ns],color='limegreen',lw=1.5,label='Experiment')
        ax1[i,j].plot(sdoft['%d' % ns], sdofe['%d' % ns],color='#FF8787', lw= 1.5, label= 'Numerical model')
        if ns == 1 or ns == 4:
            ax1[i,j].set_ylabel('Absorbed Energy (J)')
        if ns == 5 or ns == 6:
            ax1[i,j].set_xlabel('Time (s)')
        if ns == 1:
            ax1[i,j].legend(loc=2)

ax1[2,0].set_facecolor('#DEE2E6')
ax1[2,0].grid(color= 'black', ls= '--', lw= 0.5)
ax1[2,0].set_xticks([0,10,20,30])
ax1[2,0].set_yticks([0,100,200,300,400,500,600,700,800])
ax1[2,0].set_xlim([0,30])
ax1[2,0].set_ylim([0,800])
ax1[2,0].plot(timecmc['%d' % 7],energyarraycm['%d' % 7],color='limegreen',lw=1.5,label='Experiment')
ax1[2,0].plot(sdoft['%d' % 7], sdofe['%d' % 7],color='#FF8787', lw= 1.5, label= 'Numerical model')
ax1[2,0].set_ylabel('Absorbed Energy (J)')
ax1[2,0].set_xlabel('Time (s)')


fig1.delaxes(ax1[2,1])
fig1.delaxes(ax1[2,2])

fig1.tight_layout()   

plt.savefig('SDOF_Energy.png',dpi=300)


fig2, ax2 = plt.subplots(nr,nc,figsize=(13,8))

for i in range(nr-1):
    for j in range(nc):
        ns = int(3*i+j+1)
        ax2[i,j].set_facecolor('#DEE2E6')
        ax2[i,j].grid(color= 'black', ls= '--', lw= 0.5)
        ax2[i,j].set_xticks([-30,-15,0,15,30])
        ax2[i,j].set_yticks([-5000,-2500,0,2500,5000])
        ax2[i,j].set_xlim([-30,30])
        ax2[i,j].set_ylim([-5000,5000])
        ax2[i,j].plot(dispdacm['%d' % ns],accelcm['%d' % ns]*m,color='limegreen',lw=1.5,label='Experiment')
        ax2[i,j].plot(sdofd['%d' % ns]*1000, sdoff['%d' % ns],color='#FF8787', lw= 1.5, label= 'Numerical model')
        if ns == 1 or ns == 4:
            ax2[i,j].set_ylabel('Force (N)')
        if ns == 5 or ns == 6:
            ax2[i,j].set_xlabel('Displacement (mm)')
        if ns == 1:
            ax2[i,j].legend(loc=2)

ax2[2,0].grid(color= 'black', ls= '--', lw= 0.5)
ax2[2,0].set_xticks([-30,-15,0,15,30])
ax2[2,0].set_yticks([-5000,-2500,0,2500,5000])
ax2[2,0].set_xlim([-30,30])
ax2[2,0].set_ylim([-5000,5000])
ax2[2,0].set_facecolor('#DEE2E6')
ax2[2,0].plot(dispdacm['%d' % 7],accelcm['%d' % 7]*m,color='limegreen',lw=1.5,label='Experiment')
ax2[2,0].plot(sdofd['%d' % 7]*1000, sdoff['%d' % 7],color='#FF8787', lw= 1.5, label= 'Numerical model')

fig2.delaxes(ax2[2,1])
fig2.delaxes(ax2[2,2])

fig2.tight_layout()   

plt.savefig('SDOF_FroceDisp.png',dpi=300)



fig3, ax3 = plt.subplots(nr,nc,figsize=(13,8))

for i in range(nr-1):
    for j in range(nc):
        ns = int(3*i+j+1)
        ax3[i,j].set_facecolor('#DEE2E6')
        ax3[i,j].grid(color= 'black', ls= '--', lw= 0.5)
        ax3[i,j].set_xticks([0,10,20,30])
        ax3[i,j].set_yticks([-30,-15,0,15,30])
        ax3[i,j].set_xlim([0,30])
        ax3[i,j].set_ylim([-30,30])
        ax3[i,j].plot(timecmc['%d' % ns],dispdacm['%d' % ns],color='limegreen',lw=1.5,label='Experiment')
        ax3[i,j].plot(sdoft['%d' % ns],sdofd['%d' % ns]*1000,color='#FF8787', lw= 1.5, label= 'Numerical model')
        if ns == 1 or ns == 4:
            ax3[i,j].set_ylabel('Time (s)')
        if ns == 5 or ns == 6:
            ax3[i,j].set_xlabel('Displacement (mm)')
        if ns == 1:
            ax3[i,j].legend(loc=2)

ax3[2,0].set_facecolor('#DEE2E6')
ax3[2,0].grid(color= 'black', ls= '--', lw= 0.5)
ax3[2,0].set_xticks([0,10,20,30])
ax3[2,0].set_yticks([-30,-15,0,15,30])
ax3[2,0].set_xlim([0,30])
ax3[2,0].set_ylim([-30,30])
ax3[2,0].plot(timecmc['%d' % 7],dispdacm['%d' % 7],color='limegreen',lw=1.5,label='Experiment')
ax3[2,0].plot(sdoft['%d' % 7],sdofd['%d' % 7]*1000,color='#FF8787', lw= 1.5, label= 'Numerical model')

fig3.delaxes(ax3[2,1])
fig3.delaxes(ax3[2,2])

fig3.tight_layout()   

plt.savefig('SDOF_Disp.png',dpi=300)


'''
Ta = np.linspace(0,1,100)

SaE = {}
SdE = {}
SaN = {}
SdN = {}

for i in range(len(testnumber)):
    SaE['%d' % (i+1)],SdE['%d' % (i+1)] = spectrumAD(timecmc['%d' % (i+1)],accelcm['%d' % (i+1)],Ta)
    SaN['%d' % (i+1)],SdN['%d' % (i+1)] = spectrumAD(sdoft['%d' % (i+1)],sdofa['%d' % (i+1)]/1000,Ta)


fig4, ax4 = plt.subplots(nr,nc,figsize=(13,8))

for i in range(nr-1):
    for j in range(nc):
        ns = int(3*i+j+1)
        ax4[i,j].set_facecolor('#DEE2E6')
        ax4[i,j].grid(color= 'black', ls= '--', lw= 0.5)
        #ax4[i,j].set_xticks([0,10,20,30])
        #ax4[i,j].set_yticks([0,100,200,300,400,500,600,700,800])
        #ax4[i,j].set_xlim([0,30])
        #ax4[i,j].set_ylim([0,800])
        ax4[i,j].plot(SdE['%d' % ns],SaE['%d' % ns],color='limegreen',lw=1.5,label='Experiment')
        ax4[i,j].plot(SdN['%d' % ns], SaN['%d' % ns],color='#FF8787', lw= 1.5, label= 'Numerical model')
        if ns == 1 or ns == 4:
            ax4[i,j].set_ylabel('Spectral Acceleration (g)')
        if ns == 5 or ns == 6:
            ax4[i,j].set_xlabel('Spectral Displacement (m)')
        if ns == 1:
            ax4[i,j].legend(loc=2)

ax4[2,0].set_facecolor('#DEE2E6')
ax4[2,0].grid(color= 'black', ls= '--', lw= 0.5)
#ax4[2,0].set_xticks([0,10,20,30])
#ax4[2,0].set_yticks([0,100,200,300,400,500,600,700,800])
#ax4[2,0].set_xlim([0,30])
#ax4[2,0].set_ylim([0,800])
ax4[2,0].plot(SdE['%d' % 7],SaE['%d' % 7],color='limegreen',lw=1.5,label='Experiment')
ax4[2,0].plot(SdN['%d' % 7], SaN['%d' % 7],color='#FF8787', lw= 1.5, label= 'Numerical model')
ax4[2,0].set_ylabel('Spectral Acceleration (g)')
ax4[2,0].set_xlabel('Spectral Displacement (m)')


fig4.delaxes(ax4[2,1])
fig4.delaxes(ax4[2,2])

fig4.tight_layout()   

plt.savefig('SDOF_InSpectra.png',dpi=300)
'''